##Alison E. Post##
##Foreign and Domestic Investment in Argentina (2014)##
##Replication Code for Time Series Analysis (Table 7.3) and Alternative Explanations Discussion ###

library(foreign)
library(survival)
library(Zelig)
library(car)

##Load data and create relevant subset##
foreign_dom <- read.csv("replication.timeseries.csv", header=TRUE) 

##Create subset of data with just concessions and divestitures##
div.con <- foreign_dom[foreign_dom$Type_PPI!="Greenfield Project",]
div.con <- div.con[div.con$Type_PPI!="Management, Lease Contract",] 

###RESULTS FOR TABLE 7.3###
	##Note: all analyses for the book were originally conducted using the R package "Zelig".  After the proofs for the book were submitted, Zelig stopped supporting the Cox proportional hazard model. I have therefore added R survival package commands to obtain key table results throughout this file.##

##Model 1: Results similar to those from cross-sectional analysis##
z.out <- zelig(Surv(start,stop,Cancelled_In_Year) ~ Half_Foreign + GDPpc1991 + WB_Rule_Law + Polity + strata(Type_PPI), model="coxph", data=div.con, robust=TRUE, cluster="Project_name")
summary(z.out)
extractAIC(z.out)

s.out <- coxph(Surv(start, stop, Cancelled_In_Year) ~ Half_Foreign + GDPpc1991 + WB_Rule_Law + Polity + strata(Type_PPI) + cluster(Project_name), robust=TRUE, data=div.con)
summary(s.out)
extractAIC(s.out)
cox.zph(s.out)


##Model 2: Results roughly the same as cross-sectional analysis##
z.out <- zelig(Surv(start,stop,Cancelled_In_Year) ~ Half_Foreign + GDPpc1991 + WB_Rule_Law + Polity + Size_Cat + strata(Type_PPI), model="coxph", data=div.con, robust=TRUE, cluster="Project_name")
summary(z.out)
extractAIC(z.out)

s.out <- coxph(Surv(start, stop, Cancelled_In_Year) ~ Half_Foreign + GDPpc1991 + WB_Rule_Law + Polity + Size_Cat + strata(Type_PPI) + cluster(Project_name), robust=TRUE, data=div.con)
summary(s.out)
extractAIC(s.out)
cox.zph(s.out)


##Model 3: 
z.out <- zelig(Surv(start,stop,Cancelled_In_Year) ~ Half_Foreign + GDPpc1991 + Checks + FH_CL + Leftgov + Leftgov*Checks + strata(Type_PPI), model="coxph", data=div.con, robust=TRUE, cluster="Project_name")
summary(z.out)
extractAIC(z.out)

s.out <- coxph(Surv(start, stop, Cancelled_In_Year) ~ Half_Foreign + GDPpc1991 + Checks + FH_CL + Leftgov + Leftgov*Checks + strata(Type_PPI) + cluster(Project_name), robust=TRUE, data=div.con)
summary(s.out)
extractAIC(s.out)
cox.zph(s.out)


##Model 4: Results roughly the same##
z.out <- zelig(Surv(start,stop,Cancelled_In_Year) ~ Half_Foreign + GDPpc1991 + Checks + FH_CL + Leftgov + Leftgov*Checks + After.Crisis + After.Crisis*Half_Foreign + strata(Type_PPI), model="coxph", data=div.con, robust=TRUE, cluster="Project_name")
summary(z.out)
extractAIC(z.out)

s.out <- coxph(Surv(start, stop, Cancelled_In_Year) ~ Half_Foreign + GDPpc1991 + Checks + FH_CL + Leftgov + Leftgov*Checks + After.Crisis + After.Crisis*Half_Foreign + strata(Type_PPI) + cluster(Project_name), robust=TRUE, data=div.con)
summary(s.out)
extractAIC(s.out)
cox.zph(s.out)


##Interpreting the left government interaction term (p. 202, fn 32)##
	#Create data subset without missing values for simulations#
attach(div.con)
div.con.subset <- data.frame(Project_name, Lifespan, Cancelled, start, stop, Cancelled_In_Year, Half_Foreign, GDPpc1991, Country,Type_PPI, Checks, FH_CL, Leftgov, After.Crisis)
div.con.subset <- na.omit(div.con.subset)
detach(div.con)

z.out1 <- zelig(Surv(start, stop, Cancelled_In_Year) ~ Half_Foreign + GDPpc1991 + FH_CL + Leftgov + Leftgov*Checks + Checks + After.Crisis + strata(Type_PPI), model="coxph", robust=TRUE, data=div.con.subset, cluster="Project_name") 
summary(z.out1)
extractAIC(z.out1)
x.low1 <- setx(z.out1, Checks=1, strata="Type_PPI=Concession") 
s.out1 <- sim(z.out1, x = x.low1)
summary(s.out1)

x.low1 <- setx(z.out1, Checks=1, Leftgov=1, strata="Type_PPI=Concession") 
s.out1 <- sim(z.out1, x = x.low1)
summary(s.out1)


	##Graph for interpretation of interaction term coefficients##
checks.range <- 1:6
x.low <- setx(z.out1, Leftgov = 0, Checks = checks.range)
x.high <- setx(z.out1, Leftgov = 1, Checks = checks.range)
s.out <- sim(z.out1, x = x.low, x1 = x.high)
plot.ci(s.out, xlab = "Level of Checks and Balances",
ylab = "Probability of Contract Cancellation",
main = "Effect of Left government and Checks on Contract Cancellation")



###Variable Profile Table (Online appendix, table A.3)###
min(na.omit(div.con$After.Crisis))
max(na.omit(div.con$After.Crisis))
mean(na.omit(div.con$After.Crisis))
sd(na.omit(div.con$After.Crisis))

min(na.omit(div.con$Checks))
max(na.omit(div.con$Checks))
mean(na.omit(div.con$Checks))
sd(na.omit(div.con$Checks))

min(na.omit(div.con$Leftgov))
max(na.omit(div.con$Leftgov))
mean(na.omit(div.con$Leftgov))
sd(na.omit(div.con$Leftgov))

min(na.omit(div.con$FH_CL))
max(na.omit(div.con$FH_CL))
mean(na.omit(div.con$FH_CL))
sd(na.omit(div.con$FH_CL))